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Energy transport and confinement in tokamak fusion plasmas is usually determined by the cou- 
pled nonlinear interactions of small-scale drift turbulence and larger scale coherent nonlinear struc- 
tures, such as zonal flows, together with free energy sources such as temperature gradients. Zero- 
dimensional models, designed to embody plausible physical narratives for these interactions, can help 
identify the origin of enhanced energy confinement and of transitions between confinement regimes. 
A prime zero-dimensional paradigm is predator-prey or Lotka-Volterra. Here we extend a successful 
three- variable (temperature gradient; microturbulence level; one class of coherent structure) model 
in this genre [M A Malkov and P H Diamond, Phys. Plasmas 16, 012504 (2009)], by adding a 
fourth variable representing a second class of coherent structure. This requires a fourth coupled 
nonlinear ordinary differential equation. We investigate the degree of invariance of the phenomenol- 
ogy generated by the model of Malkov and Diamond, given this additional physics. We study and 
compare the long-time behaviour of the three-equation and four-equation systems, their evolution 
towards the final state, and their attractive fixed points and limit cycles. We explore the sensitivity 
of paths to attractors. It is found that, for example, an attractive fixed point of the three-equation 
system can become a limit cycle of the four-equation system. Addressing these questions which we 
together refer to as robustness for convenience, is particularly important for models which, as here, 
generate sharp transitions in the values of system variables which may replicate some key features 
of confinement transitions. Our results help establish the robustness of the zero-dimensional model 
approach to capturing observed confinement phenomenology in tokamak fusion plasmas. 

Keywords: Tokamak confinement regimes, zero-dimensional modelling, predator-prey, Lotka- 
Volterra 



1. Introduction 

Energy transport in toroidal magnetically confined fu- 
sion plasmas is determined, in most cases, by the ef- 
fects of small-scale turbulence and larger scale coherent 
nonlinear structures, together with their mutual interac- 
tions. These structures include zonal flows and geodesic 
acoustic modes [IHZ], which are radially localised poloidal 
flows, and streamers [8 , which are radially highly elon- 
gated and poloidally localised. The importance of these 
structures for energy transport was highlighted in large 
scale numerical simulations [9", TO , and the first direct 
experimental observation of streamers was reported in 
2008 [8 . Zonal flows have been the subject of extensive 
theoretical and observational work[T]-[7]. There is now 
substantial experimental support for the long-standing 
hypothesis pT] that the growth of zonal flows is driven 
by the averaged Reynolds stress of small scale turbu- 
lence. The latter can be locally suppressed by the resul- 
tant shear flow, thereby generating a temporally quasi- 
discontinuous enhancement of global energy confinement: 
the L-H transition [12]. Whether zonal flows or stream- 
ers are preferentially formed under specific plasma con- 
ditions, and how they compete, has been addressed from 
various perspectives [ 13HT5] . and remains an open ques- 
tion. For a recent review of experimental observations 
of the interaction between mesoscale structures (such 
as zonal flows and streamers) and microscale structures 
(such as drift turbulence), see [16]; of drift turbulence. 



particularly in relation to transitions in global confine- 
ment, see[17j and of the L-H transition, see[T8]. A re- 
cent review of these physics issues in a broad context is 
provided by[T9^. As emphasised in p!6UT9 | and references 
therein, recent diagnostic advances are transforming the 
experimental study of time evolving microturbulence and 
coherent nonlinear mesoscale structures during confine- 
ment transitions. This generates fresh theoretical chal- 
lenges. In addition, the ability to understand and control 
this plasma physics phenomenology will be central to the 
successful operation of the next step magnetic confine- 
ment fusion experiment ITER[20 . 

It is remarked by Malkov and Diamond in[21 , here- 
after referred to as MD, that transport models derived 
from the fundamental equations of plasma physics con- 
tinue to add much to our understanding but "tend to 
be increasingly, if not excessively, detailed. Therefore, 
there is high demand for a simple, illustrative theoreti- 
cal model with a minimal number of critical quantities 
responsible for the transition. Such models usually yield 
or encapsulate basic insight into complicated phenom- 
ena." One approach in fusion plasmas is that of zero- 
dimensional models for the interaction between micro- 
turbulence and coherent nonlinear structures, in particu- 
lar predator-prey or Lotka-Volterra [22l [23] . The proper- 
ties of Lotka-Volterra systems, both mathematically and 
from the perspective of fusion plasma physics, are by 
no means fully explored and remain an active field of 
research [24] - [29] . For fusion applications, a key step is to 
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establish agreement between the outputs of such mod- 
els and the observed confinement phenomenology, which 
should ideally extend to the character of measured time 
traces of key properties near transitions, for example. 
Recent experimental results [31, 32 are encouraging in 
this respect. There is an important additional require- 
ment. The zero-dimensional models used for this appli- 
cation should be robust, in the sense that the character 
of their outputs remains largely invariant against minor 
changes in the formulations of the models. This require- 
ment for robustness has been explicitly noted [33] in the 
other main class of zero-dimensional heuristic model for 
magnetised plasma confinement, namely sandpiles, both 
in fusion [34ti4Q] and in solar-terrestrial [33l [4TH431 J con- 
texts, and requires investigation for predator-prey and 
Lotka-Volterra applications to fusion plasmas. 

There are several aspects to the degree of invariance 
of the phenomenology generated by a zero-dimensional 
model when aspects of the model are changed. First, 
what is the long-time behaviour of the system and how 
sensitive is this to variation in the model parameters [44", 
[45j? Second, how sensitively does the nature of the sys- 
tem's evolution towards its final state depend on the ini- 
tial conditions? Is there an attractive fixed point or limit 
cycle towards which the system fiows as time passes? If 
so, what is its basin of attraction? Third, how sensitive is 
the path to this attractor? This is particularly important 
for models which, as here, generate sharp transitions in 
the values of system variables which may replicate some 
key features of confinement transitions in tokamaks. If 
the initial conditions are varied, is the time at which the 
transition occurs delayed or brought forward, or does its 
character change, for example? Further, given two zero- 
dimensional models which are schematically distinct but 
adjacent, how similar is the phenomenology of their solu- 
tions? An example is provided here by our extension of 
the model of MD[21 to incorporate two variables, rather 
than one, representing different classes of large scale co- 
herent nonlinear field, in a four- variable system. The case 
of two predators and one prey was considered theoreti- 
cally in the model of Itoh & Itoh^29j, hereafter referred 
to as II, and by Miki and Diamond [30], and there is re- 
cent experimental motivation [3T| 132], Insofar as a zero- 
dimensional model turns out to be robust with respect 
to the considerations outlined (attractors; initial condi- 
tions; structural adjacency), confidence is strengthened 
in the mapping from model variables to specific plasma 
properties, and from the time evolving behaviour of the 
model to that of the plasma system. 

In the present paper, we focus from this perspective on 
the interesting and successful mathematical model pro- 
posed in MD. This is constructed in terms of variables 
representing the magnitude of the plasma temperature 
gradient and the amplitudes of small scale drift turbu- 
lence and of large scale coherent nonlinear structures 
such as zonal flows. Malkov & Diamond proposed [21] 
certain mappings between different solution regimes of 
their model and different confinement regimes of toka- 



mak plasmas. In the interest of continuity, we follow the 
confinement regime nomenclature of MD in relation to 
model outputs in the present paper. We investigate the 
robustness of the phenomenology of the MD model ex- 
tended as described, for parameter regimes identical, or 
adjacent, to those used in the key figures of MD. Where 
robustness is demonstrated and, if possible, explained, 
this reinforces confidence that models in the genre of MD 
and II may capture key features of the physics of confine- 
ment transitions in tokamak plasmas. 

2. Model equations 

Specifically, the MD model is a closed system of nonlin- 
ear differential equations which couple the time evolution 
of three variables: the drift wave-driving temperature 
gradient TV, the energy density of drift wave turbulence 
and the zonal fiow velocity U. The three variables 
of the II model exclude TV, retain drift turbulence en- 
ergy density denoted by and incorporate the energy 
densities of two competing classes of coherent nonlinear 
structure, zonal fiows Z and zonal fields (e.g. streamers) 
M. Miki and Diamond [30 introduced a zero-dimensional 
three-variable two-predator, one prey model, where the 
predators are identified with zonal fiows and geodesic 
acoustic modes. The aspect of robustness which we first 
address can therefore be expressed in physical terms as 
follows. We adopt the philosophy of II and of Ref.[30 by 
introducing two competing classes of coherent nonlinear 
structure, here identified with zonal flows and streamers, 
that replace the single class in MD. The other two MD 
equations are adjusted only so far as necessary to accom- 
modate these two flelds, instead of one, in a mathemati- 
cally symmetrical way as in II. We investigate how far the 
model outputs of our new four- variable system differ from 
those of the three- variable system of MD. A good focus 
for this study is provided by the time traces captured 
in Figs. 2-4 of MD, which have been mapped to tran- 
sitions observed between tokamak conflnement regimes. 
How are these traces altered by the inclusion of a second 
competing class of coherent nonlinear structure? The 
answers to these questions are conditioned by the un- 
derlying phase space structure of families of solutions to 
the models, as plotted in Fig. 5 of MD, for example. In 
addition to studying time traces, therefore, we seek to 
characterise the limit cycles and flxed points of our sys- 
tem of equations. We flrst generalize the un- normalized 
MD equations to: 
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This model encompasses drift wave turbulence level 5, 
drift wave driving temperature gradient A/", zonal flow ve- 
locity VzF^ streamer flow velocity Vsf^ and the heating 
rate q which is a control parameter of the system. This 
model thus extends, to the case when zonal flows are 
joined by streamers, the key physics encapsulated in the 
description in[46 : "When the drift wave turbulence drive 
becomes sufficiently strong to overcome flow damping, it 
generates zonal flows by Reynolds stress. Drift wave tur- 
bulence and zonal flows then form a self-regulating sys- 
tem as the shearing by zonal flows damps the drift wave 
turbulence." We note that this model follows the ap- 
proach expressed in Eq.(17) of MD[21 , in that the zonal 
flows and streamers do not explicitly enter the time evolu- 
tion equation for the temperature gradient, Eq.(4). The 
zonal flows and streamers are indirectly coupled to each 
other through the evolving temperature gradient and mi- 
croturbulence level. To maximise mathematical congru- 
ence with the model of MD, there is no direct cross term 
in VsfVzf- 

The corresponding normalized equations are 
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We note that this is equivalent to the normalization of 
MD only when a2 = 1. This reflects a minor inconsis- 
tency in the normalization choice of MD. The system of 
Eqs.(5-8) thus generalizes the system of Eqs. (15-17) of 
MD by introducing two distinct flow variables, Ui and 
C/2, to replace the single zonal flow variable U. We refer 
to /7i as zonal flow, /72 as streamer flow. 

Section 3 of this paper addresses transition phe- 
nomenology given time-independent coefficients, as char- 
acterised primarily by time traces. This requires care- 
ful comparison with the specific scenarios identified in 
Fig. 3 to Fig. 5 of MD. The MD scenarios predetermine the 
choice of parameter values and initial conditions that we 
consider. We typically probe neighbouring phase space 



by considering in addition eighty-one (three to the fourth 
power) nearby phase trajectories. In Section 4 we con- 
sider the phase space evolution of our system and es- 
tablish comparisons between the MD model and ours. In 
Section 5 we analyse possible links to the phenomenology 
of tokamak plasmas, in the spirit of MD and II. 

3. Modelling confinement transitions 

In the limit where either one of the two parameters that 
represent distinct classes of coherent nonlinear structures 
(zonal flows or streamers) in our model vanishes, it re- 
produces exactly the results shown in Fig. 2 of MD, as 
required. Figure 1 displays the corresponding results for 
the case where both streamers and zonal flows exist. In 
the nomenclature of MD, the system starts from an over- 
powered state near H-mode, with negligible turbulence E 
and large scale structures [/i, U2' The eventual growth 
of turbulence accompanies a sharp drop in N to unsta- 
ble L-mode, while also providing energy for Ui and U2- 
Drift wave turbulence is later suppressed and the maxi- 
mum amplitude of large scale flows declines, leaving only 
the mean flow to support the transport barrier [19]. Fi- 
nally the stable T-mode, which combines a steady- state 
level of E with lower N than H-mode, appears after the 
oscillating transition regime. During this transition, en- 
ergy is extracted from the initially dominant oscillating 
streamer flow U2 the zonal flow Ui until the former 
vanishes. 

In Fig. 2, we plot the system evolution for the case 
where the values of V2 and 772 are different from Fig.l, 
while all other parameter values are identical. Specifi- 
cally, in Fig.l z^2/^i = ^2/^1 = 1-01, whereas in Fig. 2 
^2/^1 = 0.01 and 7^2/^1 = 0.1. This weakens both 
the drive and the damping of structures U2 compared 
to zonal flows /7i in Fig. 2, with respect to the case of 
Fig.l. Before time reaches t ~ 6000, the evolution is 
very similar to Fig. 2 of MD. However, at t~ 6500 we 
find a dramatic change. A limit cycle appears after the 
long-term fixed point time series. The amplitudes oiUi 
and U2 exchange rather fast compared to Fig.l. Further- 
more, the period of the limit cycle is rather long: several 
hundred time units. With the appearance of zonal flows 
and streamers, the T-mode becomes unstable. 

Figure 3 shows the case where the heating rate is higher 
than for Fig.l, q = 0.58, but all other model parameters 
are the same. At each pulsed occurrence of zonal flows 
Ui and streamers 1/2^ the former extract energy from 
the latter, which become extinct after the sixth pulse. 
Thereafter there are limit cycle oscillations in E^ N and 
Ui equivalent to the limit cycle for E^ N and U in the 
case in MD. 

Figure 4 shows time traces for the case where all pa- 
rameters, except the heating rate q = 0.58 which is the 
same as in Fig. 3, are those of Fig. 2. Together with Fig. 5, 
where the heating rate q is slightly increased to q = 0.582 
instead of q = 0.58, this enables us to relate our model 
to Fig.4 of MD, which showed that if in MD q = 0.582 
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Figure 1. Upper panel: From overpowered H-mode to unsta- 
ble L-mode then to T-mode. Lower panel: Transition to T- 
mode for Ui and U2 showing intersection at t 750 followed 
by energy reversal. The parameters are ui = 19, 1^2 = l.Olz/i, 
ryi = 0.12, r]2 = l.Olryi, q = 0.47, p = 0.55, a = 0.6, C = 1.7. 



Figure 2. Upper panel: Transition from stable fixed point 
state to unstable oscillatory limit cycle state. Lower panel: 
Zoom in version from t = 300 to t = 800. The parameters 
are ui 19, 1^2 O.Olz/i, r]i 0.12, r]2 O.lryi, q 0.47, 
p = 0.55, cr = 0.6, C = 1.7. 



instead of 0.58, the limit cycle eventually collapses after 
many oscillations. The final state has N finite and the 
remaining variables are zero; this is designated the QH- 
mode fixed point in MD. The corresponding cases for our 
model Eqs.(5-8) are shown in Figs. 4 and 5. A precursor 
to limit cycle collapse is apparent in Fig. 4 in the growth 
of the streamer field U2 during the episodes of zonal flow 
quiescence in the last few oscillations of the system. 

For the slightly different parameter set used to gener- 
ate Fig. 5, the pulses of /7i and U2 grow and die together. 
Their peak amplitude increases at each successive cycle, 
as does the time interval between them. At the final oscil- 
lation, Ui and U2 collapse promptly together, whereas E 
survives longer until it is extinguished by damping. The 
phenomenology of Fig. 5 thus corresponds more closely to 
that of Fig. 4 of MD, compared to our Fig. 4. 

Figure 6 illustrates how system evolution towards the 
finite- final state of Fig. 5 depends on the damping rate 
r]2 of streamers. We fix all parameters except r]2 and 
find that, with increasing 772, there are more peaks of 
U2 correlating with cyclic growth of which acts as a 
damping sink of N. Successive peaks increase in height 
prior to extinction, which results in a final state similar 
to Fig. 5. 

4. Phase space evolution 

The time traces of the individual variables, plotted in 
Figs.l to 6, represent projections of the evolution in four- 
dimensional phase space of the system defined by Eqs.(5) 
to (8). In the present section, we capture the global phase 
space explored by this system, for parameter values cor- 
responding, or adjacent, to those used to generate Figs.l 
to 6. This approach enables us to identify and charac- 
terise the nature of initial and final states, and of the 
transitional behaviour between them. These results are 
supplemented in the Appendix by stability studies. At 
issue are two main physical concerns, which map directly 
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Figure 3. Energy transfer from IJ2 to IJ\ during pulses of 
strong nonlinear oscillation, followed by limit cycle ocillation 
in A", E and U\. The parameters are v\ — 19, V2 — l.Olz^i, 
ryi = 0.12, T]2 = l.Olryi, q = 0.58, p = 0.55, a = 0.6, C = 1-7. 
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Figure 4. Upper panel: Collapse of limit cycle in A, £^ and 
IJ\. Lower panel: Stair increasing of IJ2- The parameters 
are 19, V2 O.Olz/i, ryi 0.12, 7/2 O.Olryi, q 0.58, 

p = 0.55, cr = 0.6, C = 1.7. 
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Figure 5. Upper panel: Collapse of limit cycle with positively 
correlated growth of pulses of Ui and 1/2- Lower panel: Zoom 
in version from t = 240 to t = 400. The parameters are 
z/i = 19, = l.OOOlz/i, r]i = 0.12, r]2 = l.OOOlryi, q = 0.582, 
p = 0.55, cr = 0.6, C = 1.7. 




Figure 6. Evolution to the finite N attractor for different 
values of 772. Upper panel: 772 = 0.05. Middle upper panel: 
772 = 0.06. Middle lower panel: 772 = 0.10. Lower panel: 
7^2 = 0.11. The remaining parameters are the same: ui = 19, 
ly^ = l.OOlz/i, ryi = 0.12, q = 0.582, p = 0.55, a = 0.6, C = 1-7. 



to the properties of different energy confinement regimes 
in tokamaks, insofar as the zero-dimensional approach 
and the identifications made in MD, for example, may 
be valid. First, what is the nature of the final state that 
is reached at long times? For example, is it an attractive 
fixed point or a limit cycle (implying a nearby repulsive 
fixed point)? Second, there is the question, discussed 
previously, of robustness of three- variable models against 
the inclusion of a fourth variable (here, streamers) in the 
model. For example, the pioneering work of MD includes 
identification of a limit cycle (Fig. 3 of MD) with a specific 
confinement regime. Is this limit cycle - and, proceeding 
by analogy, the confinement regime that it represents - 
stable against the presence of streamers in addition to 
zonal flows? 

Figure 7 displays the generalisation, to the four- 
variable system, of the case of the three-variable system 
addressed in Fig. 2 of MD. To fix ideas, the two left-hand 
plots correspond to the three-variable case for the pa- 



rameters of Fig. 2 of MD, showing the attractive fixed 
point which has finite values of N and U. The inward 
spiral path of the system from a random initial position 
is shown, both in {E, U) space and projected onto 
the (£^, U) plane. It is evident that this path lies on a 
topological structure in phase space, whose dimension- 
ality is lower by one than that of the full phase space. 
The two right-hand plots of Fig. 7 show how this system 
changes when the two variables Ui and U2 replace [/, for 
the parameter values used to generate the traces in Fig.l, 
which are adjacent to those for Fig. 2 of MD, as discussed 
above. The centre right-hand plot shows initial spiral 
convergence in (^, U2) which closely resembles that in 
the {E^ U) plane displayed at centre left. Whereas with 
three variables this convergence is towards a fixed point, 
the existence of a fourth variable renders this attractive 
fixed point unstable. In consequence, the final stage of 
system evolution consists of injection in the Ui direc- 
tion to a fixed point at finite {E^ N ^ Ui) with U2 = 0. 
The far right plot in Fig. 7 demonstrates that this is in- 
deed a fixed point, towards which phase space evolution 
originating from eighty-one different initial points con- 
verges. In each case, there is spiral convergence on a 
manifold followed by injection along Ui. The choice of 
initial condition affects only the orientation of this con- 
vergence manifold with respect to Ui and U2. We note 
also that the final state with finite Ui differs from the 
MD final state for which U = 0. 

Figure 8 illustrates the phase space evolution of the 
system whose time traces are plotted in Fig. 2, which like 
Fig. 7 is a case with parameters adjacent to those used to 
generate Fig. 2 of MD. The initial spiral convergence in 
the {E^ Ui) plane, shown in the centre panel, resembles 
that in the {E^ U) plane for the MD case plotted in the 
left panel, which is identical to the centre-left panel of 
Fig. 7. As in Fig. 7, the stable fixed point of the three- 
variable system is unstable for the four-variable system, 
for which there is injection along 112- Unlike Fig. 7, where 
this injection is towards a stable fixed point, in Fig. 8 the 
injection is onto a stable limit cycle that has finite slow 
oscillations in (A/", E^ U2) with Ui = in the four- variable 
system. 

The three- variable MD system has a limit cycle in (A', 
E^ U) for the case shown in Fig. 3 of MD. This is re- 
plotted in the two left panels of Fig. 9 and in the left 
panel of Fig. 10. Figures 9 and 10 relate to the time traces 
shown in Figs. 3 to 5 of the present paper, obtained for 
parameter sets for the four- variable system which are ad- 
jacent to those used in MD for the three- variable system. 
For the parameters of Fig. 9, which is the phase space plot 
for Fig. 3, it is clear from the two right-hand panels that 
the limit cycle behaviour is essentially that of the MD 
system. The transient evolution towards the limit cycle 
involves circulation on similar planes that have succes- 
sively lower peak values oi U2- The final limit cycle in 
(A', E,Ui), with U2 0, is essentially that in (A', E, U) 
for the three- variable system. 

The three-variable MD attractive limit cycle which 
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manifests in the four- variable system as shown in Fig. 9 
is, however, unstable. Figure 10, which is the phase space 
plot for Fig. 4, shows that the system leaves the former 
limit cycle and transiently explores the additional phase 
space dimension associated with the additional variable, 
before converging to a new fixed point that has N finite 
and all other variables zero. This class of attractive fixed 
point is noted in Fig. 4 of MD, shown in the far left panel 
of Fig. 11 and, projected on the {E^ U) plane, in the cen- 
tre left panel. The two right-hand panels of Fig. 11 are 
the phase space plots for Fig. 5, showing convergence to 
the origin in /7i, U2) space while N remains finite. 
The final step to the origin is preceded by circulation 
around and away from an apparent repulsive fixed point 
with finite values of Ui and The far right panel of 
Fig. 11 shows that the choice of initial conditions merely 
affects the orientation in ([/i, U2) space of the plane of 
this transient circulation. 

The phase space behaviour discussed thus far assists 
us in re- visiting the time traces in Fig. 2, for which the 
corresponding phase plot is given in Fig. 13. In Fig. 12 
we annotate Fig. 2 in light of Fig. 13. These two Figures 
demonstrate how, for the four-variable system, the T- 
mode of the three-variable system becomes unstable at 
long times. The system then evolves towards the newly 
identified attractive limit cycle in (A/", 1/2)- Here slow 
oscillations in N correlate with those in 1/2^ both of which 
remain finite throughout, while bursts of E, feeding on 
[/2, occur between extinctions. 

5. Conclusions 

Contemporary experimental results from the DIII- 
D[3T| and HL-2A tokamaks[32 reinforce the relevance of 
zero-dimensional predator-prey models to transitions be- 
tween energy confinement regimes. Understanding how 
the outputs of related, but different, predator-prey mod- 
els for plasma confinement phenomenology may resemble 
or deviate from each other is therefore important. In this 
paper we have focused on the consequences of adding a 
second predator, and hence a fourth field variable, to 
the three-field MD[21 model. Quantitative studies have 
been presented for parameter sets that are maximally ad- 
jacent to those in MD, which yield the time traces shown 
in Figs.l to 6 and Fig. 12. These are projections of the 
phase space dynamics shown in Figs. 7 to 11 and Fig. 13. 
It is found that both congruences and deviations can oc- 
cur between the three-field and four-field models. For ex- 
ample. Fig. 10 shows how a limit cycle in the three-field 
system is unstable for four fields in the relevant parame- 
ter range, where the attr actor is a fixed point. Conversely 
Fig. 8 shows a three-field fixed point mapping to a four- 
field limit cycle. Figure 13 shows the complex, but re- 
solved, phase space dynamics underlying a generalisation 
to four fields of the three-field scenario modelled in Fig. 2 
of MD. We conclude that exploration of the linkages be- 
tween different zero-dimensional models, capturing full 
phase space properties so far as computationally possi- 



ble, needs to keep pace with the continuing development 
and refinement of individual zero-dimensional models in 
fusion plasma physics. 

Zero-dimensional models remain attractive because 
they embody physically motivated narratives that may 
account for global fusion plasma confinement phe- 
nomenology. Ideally the end states (attractors) of zero- 
dimensional models, together with the transitional be- 
haviour en route from the initial configurations, should 
be robustly identifiable with fusion plasma confinement 
states and transitions. Zero-dimensional predator-prey 
models, constructed in terms of a small number of vari- 
ables representing global quantities such as the drift wave 
turbulence level f , drift wave driving temperature gra- 
dient A/", zonal fiow velocity Vzf^ streamer fiow velocity 
VsF^ and the heating rate q in Eqs.(l) to (4), are intrin- 
sically nonlinear. This nonlinearity implies the potential 
for a rich and varied set of attractors and transitional 
behaviour, together with strong dependence on the nu- 
merical values of model parameters. The present paper 
has taken steps to explore this potential for the model of 
interest in the case of parameter sets close to those stud- 
ied previously in MD, with a view to strengthening the 
links between families of zero-dimensional models on the 
one hand, and fusion plasma confinement phenomenology 
on the other. We note finally that some of the considera- 
tions addressed here may carry over to other fields where 
it is hoped to develop zero-dimensional models that have 
descriptive, or even predictive, power for global phenom- 
ena in macroscopic multiscale driven-dissipative systems. 
A topical instance is provided by zero-dimensional mod- 
elling in climate science, see for example Ref.[47 and ref- 
erences therein, where some general circulation models 
incorporate Lotka-Volterra features [48]. 
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APPENDIX: Identification and stability of fixed 
points 

We start from Eqs.(5-8), and for simplicity define the 
normalized equations as 

'dE/dt = {N - - E - Ui - U2) E = f {E, Ui, t/2, N) 
dUi/dt = i/i r - vi) Ui = gi {E, Ui,N) 

dU2/dt = z/2 r - m) U2 = 92 {E, C/2, N) 

JN/dt = q-lp^aE)N = h {E, N) 

(9) 

We regard point (A^o^ ^o? t/10, U20) as a fixed point of 
the 4D system, and define 

' fo = f{Eo,Uw,U2o,No) 
^ 910 = 9i(Eo,Uio,No) ^^Q^ 
520 = 92 {Eo,U2o,No) 
^ho = h{Eo,No) 
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Figure 7. First panel: Fig.2 in MD. The parameters are = 19, 77 = 0.12, q = 0.47, p = 0.55, a = 0.6, C = 1-7. Second 
panel: Projection of first panel on E-U plane. Third panel: Phase plot of Fig.l. Last panel: Phase plot of Fig.l with 81 initial 
conditions. Stars denote initial values, blue dots denote trajectories and red diamonds denote final states. 




Figure 8. First panel: Projection of Fig.2 in MD on E-U plane. The parameters are = 19, 77 = 0.12, q = 0.47, p = 0.55, 
a — 0.6, — 1.7 . Second panel: Phase plot of Fig.2. Last panel: Phase plot of Fig.2 with 81 initial conditions. Stars denote 
initial values, blue dots denote trajectories and red diamonds denote final states. 
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Figure 9. First panel: Fig. 3 in MD. The parameters are = 19, 77 = 0.12, q — 0.58, p — 0.55, a — 0.6, C — 1-7. Second 
panel: Projection of first panel on E-U plane. Third panel: Phase plot of Fig. 3. Last panel: Phase plot of Fig. 3 with 81 initial 
conditions. Stars denote initial values, blue dots denote trajectories and red diamonds denote final states. 




Figure 10. First panel: Projection of Fig. 3 in MD on E-U plane. The parameters are = 19, 77 = 0.12, q = 0.58, p = 0.55, 
a — 0.6, C — 1.7. Middle panel: Phase plot of Fig. 4. Last panel: Phase plot of Fig. 4 with 81 initial conditions. Stars denote 
initial values, blue dots denote trajectories and red diamonds denote final states. 
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Figure 11. First panel: Phase plot for Fig. 4 of MD. Second panel: Projection of Fig. 4 in MD on E-U plane. The parameters 
are = 19, ?7 = 0.12, q = 0.582, p = 0.55, a = 0.6, ( = 1.7. Third panel: Phase plot of Fig. 5. Last panel: Phase plot of Fig. 5 
with 81 initial conditions. Stars denote initial values, blue dots denote trajectories and red diamonds denote final states. 



Case 


q 




V2/V1 


Timetraces 


Phaseplot 


Manifold 


1 


0.47 


1.01 


1.01 


Fig.l 


Fig.7 


Fixed point 


2 


0.47 


0.01 


0.1 


Fig.2 


Fig.8 


Limit cycle 


3 


0.58 


1.01 


1.01 


Fig.3 


Fig.9 


Limit cycle 


4 


0.58 


0.01 


0.01 


Fig.4 


Fig. 10 


Limit cycle 


5 


0.582 


1.0001 


1.0001 


Fig.5 


Fig. 11 


Fixed point 


6 


0.582 


1.001 


0.05;0.06;0.1;0.11 


Fig.6 


N/A 


N/A 



Table L Summary of Figs.l to 11 



By construction /o = gio = ^20 = ho = 0. Near the 
fixed point, we make a local linear expansion of the model 
parameters: 



AE = E-Eo; AUi = Ui-Uio] A/72 = U2-U20; 
This gives rise to the linearized equations 



AN = N- 
(11) 



-Nn 



92 



/o4 

- 920 



^ AE ^ 
dE ^ ^ ^ 

^ dE ^ ^ 



AUi 



df 
dUi 



^ AN 



(12) 

To obtain the eigenvalues of the system, we calculate 
the corresponding Jacobian matrix 



dEJ 

9 n 

dE^l 
dE92 



d 
dUi 

dlhd2 



f dUo f 



\ dE'^ dUi"" dUo 



d 

dU2^ 

h 



dU2' 



m92 



(13) 



{Eo,Uio,U20,Nq) 



We now identify the fixed points. 
(T)ifE = 0, 



N ~N^■ 
Ul=0 

U2=0 



E — Ui — U2 = Anyvalue 



(14) 



11 



(N -N^-E-Ui-U2 = 
E 

^ - 7/1 ) [/i = 



< 



1 + C^4 

E 



U-{p^crE)N = ^ 



(15) 



From the second and third equations in this group, it 
fohows that Ui and U2 cannot be non-zero simultane- 
ously. 

(i) if Ui=0,U2^Q,E^ 0, 



Solutions for the specific cases of the MD and ZCD 
systems considered in this paper are shown in Tables II 
and Table III. 



0^ [lOOO 



2000 3000 4000 5000 



9000 10000 



situation 1 

situation 4 




N 

— E 
UI 

— U2 







300 350 400 450 500 550 



700 750 



N - - E -Ui-U2 = ^ 

^ _^ - ^1 = Anyvalue 
E 



l + CA^4-^2 = 
^q-{p^(jE)N = {) 



(16) 



Figure 12. Time series of Fig. 2 in this paper, annotated in 
light of Fig. 13. 



(ii) if /7i 7^ 0, [72 = 0, ^ 0, 



N - - E -Ui-U2 = ^ 
E 

Vi=0 



7V4 



r]2 — Anyvalue 



^q- {p^aE)N = 



(iii) ifUi = U2=0,E^ 0, 



(17) 



Qg-_ "-^ fixed point 3 
0.7- I ' '■ ' — 



Figure 13. Phase plot of Fig. 2 in this paper. 



'N-N*-E = 
Ui=0 

U2=0 

^q-{p + aE)N = 



(18) 



12 



MD 


Fixed points 


Eigenvalues 


Property 


Fig.2 


E = 0',U = 0',N = 0.8545 


-2.28;-0.55;0.3213 


Saddle point-Index 1 


E = 0.1742;t/ = 0.2780;A^ = 0.7181 


-0.0360 zb 0.8099z;-0.7567 


Spiral node (final state) 


E = 0.4638;[/ = 0;A^ = 0.5675 


-0.6460 zb0.0963z;5.2111 


Inward spiral and source 


Fig.3 


E = 0;U = 0;N = 0.8545 


-2.28;-0.55;-0.1821 


Node 


E = 0.2249;t/ = 0.1077;iV = 0.8468 


0.0069 zb0.499H;-0.9236 


Outward spiral and sink (limit cycle) 


E = 0.0769;C/ = 0;A^ = 0.9729 


0.0969;-0.7700;-1.7010 


Saddle point-Index 1 


E = 0.4588;[/ = 0;iV = 0.7028 


3.8817;-0.3122;-0.9718 


Saddle point-Index 1 


Fig.4 


E = 0;U = 0]N = 1.0582 


-2.28;-0.55;-0.1957 


Node (final state) 


E = 0.2260;t/ = 0.1036;A^ = 0.8489 


0.0080 zb 0.4892z;-0.9275 


Outward spiral and sink 


E = 0.0825;t/ = 0;iV = 0.9708 


0.1002;-0.7821;-1.6558 


Saddle point-Index 1 


E = 0.4576;C/ = 0]N = 0.7058 


3.8348;-0.3058;-0.9764 


Saddle point-Index 1 



Table II. MD system 



ZCD 


Fixed points 


Eigenvalues 


Property 


Fig.7 


E = 0;Ui = 0;U2 = 0;N = 0.8545 


-2.28;-0.55;0.3213;-2.3258 


4D Saddle point-Index 1 


E = 0.1757;t/i = 0]U2 = 0.2770;7V = 0.7171 


-0.0365 zb 0.8165z;-0.7581;0.0228 


Inward spiral, source and sink 


E = 0.1742;C/i = 0.2780;C/2 = 0;A^ = 0.7187 


-0.0360 zb 0.8099z;-0.0230;-0.7567 


Spiral node (final state) 


E = 0.4638;C/i = 0;U2 = 0;N = 0.5675 


-0.6460 zb 0.0963z;5.2402;5.2111 


Inward spiral and sources 


Fig.8 


E = 0;[/i = 0;U2 = 0;N = 0.8545 


-2.28;-0.55;0.3213;-0.0023 


4D Saddle point-Index 1 


E = 0.0219;t/i = 0]U2 = 0.3275;A^ = 0.8346 


0.0019 zb 0.0272z;-0.5888;-2.052 


Outward spiral and sinks (limit cycle) 


E = 0.1742;C/i = 0.2780;C/2 = 0;A^ = 0.7187 


-0.0360 zb 0.8099z;0.0205;-0.7567 


Inward spiral, source and sink 


E = 0.4638;t/i = 0;U2 = 0;A^ = 0.5675 


-0.6460 zb 0.0963i;0.0726;5.2111 


Outward spiral and sources 


Fig.9 


E = 0;C/i 0;U2 = 0;iV = 1.0545 


-2.28;-0.55;-2.3258;-0.1821 


Node 


E = 0.2265;[/i = 0;U2 = 0.1078;7V = 0.8456 


0.0062 zb 0.5045z;0.0228;-0.9248 


Outward spiral, source and sink 


E = 0.2249;t/i = 0.1077;t/2 = 0;A^ = 0.8468 


0.0069 zb 0.4991z;-0.0230;-0.9236 


Outward spiral and sinks (limit cycle) 


E = 0.0769;t/i = 0;U2 = 0;iV = 0.9729 


0.0969;-0.7700;-1.7010;-1.741 


4D Saddle point-Index 1 


E = 0.4588;t/i = 0;t/2 = 0;N = 0.7028 


-0.3122;-0.9718;3.8817;3.8975 


4D Saddle point-Index 2 


Fig. 10 


= 0;Ui = 0;U2 = 0;N = 1.0545 


-2.28;-0.55;-0.1821;-0.0002 


Node (final state) 


E = 0.2249;t/i = 0.1077;t/2 = 0;iV = 0.8468 


0.0069 ± 0.4990i;0.0226;-0.9236 


Outward spiral, source and sink 


E = 0.4588;t/i = 0;U2 = 0;iV = 0.7028 


-0.3122;-0.9718;3.8817;0.0614 


4D Saddle point-Index 2 


E = 0.0769;t/i = 0;t/2 = 0;iV = 0.9729 


0.0969;0.0056;-0.7700;-1.7010 


4D Saddle point-Index 2 


Fig. 11 


E = 0;C/i = 0;C/2 = 0]N = 1.0582 


-2.28;-0.55;-2.2805;-0.1957 


Node (final state) 


£^ = 0.2260;[/i = 0;[/2 = 0.1036;iV = 0.8489 


0.0080 zb 0.4892i;0.0002;-0.9275 


Outward spiral, source and sink 


E = 0.2260;t/i = 0.1036;t/2 = 0;A^ = 0.8489 


0.0080 zb 0.4892z;-0.0002;-0.9275 


Outward spiral and sinks 


E = 0.0825;C/i = 0;C/2 = 0;A^ = 0.9708 


0.1002;-0.7821;-1.6558;-1.6562 


4D Saddle point-Index 1 


E = 0.4576;C/i = 0;C/2 = 0;iV = 0.7058 


-0.3058;-0.9764;3.8348;3.8349 


4D Saddle point-Index 2 



Table III. ZCD system 



